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Astrophysical accretion is arguably the most prevalent physical process in the Universe; it occurs 
during the birth and death of individual stars and plays a pivotal role in the evolution of entire galaxies. 
Accretion onto a black hole, in particular, is also the most efficient mechanism known in nature, con- 
verting up to 40% of accreting rest mass energy into spectacular forms such as high-energy (X-ray and 
gamma-ray) emission and relativistic jets. Whilst magnetic fields are thought to be ultimately respon- 
sible for these phenomena, our understanding of the microphysics of MHD turbulence in accretion flows 
as well as large-scale MHD outflows remains far from complete. We present a new theoretical model for 
astrophysical disk accretion which considers enhanced vertical transport of momentum and energy by 
MHD winds and jets, as well as transport resulting from MHD turbulence. We also describe new global, 
3D simulations that we are currently developing to investigate the extent to which non-ideal MHD 
effects may explain how small-scale, turbulent fields (generated by the magnetorotational instability - 
MRI) might evolve into large-scale, ordered fields that produce a magnetized corona and/or jets where 
the highest energy phenomena necessarily originate. 
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1. Introduction 

It is widely accepted that high-energy astrophys- 
ical sources such as active galactic nuclei (AGN), 
gamma-ray bursts and some X-ray binaries are pow- 
ered by accretion of matter onto a central black hole. 
Since the standard theory of astrophysical disk ac- 
cretion was formulated over 30 years ago [TJ [5], ar- 
guably the most important advance in our under- 
standing of the process by which matter in the disk 
can shed its angular momentum and release its grav- 
itational binding energy has come from computa- 
tional modelling. Numerical simulations demonstrate 
unequivocally that the magnetorotational instability 
(MRI, [21 HI O E] ) can produce magnetohydrodynamic 
(MHD) turbulence and enhanced angular momentum 
transport (see [7] for a review) . The presence of even 
a very weak magnetic field is the key ingredient: it 
completely changes the dynamics from a keplerian 
flow which is hydrodynamically stable even at high 
Reynolds numbers (as recently verified experimentally 
[5]) to one which is unstable to the rapid growth of 
MHD modes leading to turbulence in the nonlinear 
regime. 

It is over 20 years since the first MHD simulations 
of astrophysical accretion flows were carried out [§1 
[TU]. Notwithstanding the important advances made 
to date [3 HH E21 H3 HH ESI Q51 113 ESI E2I , numerical 
simulations have so far been unable to resolve two key 
outstanding issues: 
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1. How are the high rates of mass accretion inferred 
in the most powerful sources achieved? 

2. How are the outflows and jets observed across the 
mass spectrum of accreting sources produced? 

In what follows, we briefly address each of these open 
questions and suggest how they may be connected and 
mutually resolved by a generalized model for MHD 
disk accretion. 

The most powerful accreting sources (i.e. quasars 
and other active galaxies) are fuelled by accretion 
onto a supermassive (1O 6_9 M ) black hole. They 
produce radiative luminosities that can exceed those 
of normal galaxies by several orders of magnitude 
(e.g. up to 10 48 ergs _1 ), indicating mass accretion 
rates which can exceed 100 solar masses per year 
(lM Q yr _1 w 6 x 10 25 gs _1 ). This is strictly a lower 
limit because accretion can also drive mechanical out- 
flows, in some cases with inferred kinetic powers that 
are considerably greater than the observed radiative 
luminosity (e.g. the famous M87 jet - see [5D] and ref- 
erences therein). Indeed, the fact that collimated jets 
are observed across a wide range of accreting sources 
(see [21] for a review), including those that are non- 
relativistic (e.g. protoplanetary systems, young stel- 
lar objects, see Fig. [T] and neutron stars, see espe- 
cially [H]), suggests that accretion provides an effec- 
tive, generic mechanism for powering jets and other 
energetic plasma outflow phenomena. 

Magnetic fields are required to collimatc jet 
plasma and to account for the observed radio syn- 




Fig. 1 Hubble Space Telescope optical images of jets in 
young stellar objects. Photo credits: C. Burrows, 
J. Morse (STScI), J. Hester (AZ State U.), NASA. 



chrotron jet emission. MHD simulations have revealed 
that a net poloidal component of the magnetic field, 
B z , is required to produce winds and jets from an 
accretion disk (see [23] fc> r a review). Interestingly, 
high mass accretion rates can only be achieved in the 
simulations if a large-scale MHD outflow is present 
(e.g. [121 EI])- It is particularly noteworthy that 
comparably high accretion rates cannot be achieved 
through MRI-driven MHD turbulence alone, despite 
the fact that the effective turbulent viscosity is "en- 
hanced" relative to kinematic fluid shear viscosity (e.g. 
[131 El [H]). Therefore, magnetized jets (and by im- 
plication, MHD outflows in general) must be primar- 
ily responsible for truly enhanced transport in astro- 
physical accretion flows. This conclusion is consistent 
with early analytical models [2H Hi] proposing that 
accretion is facilitated by the vertical transport of an- 
gular momentum resulting from MHD torques on the 
disk surface. It is also consistent with non-ideal MHD 
simulations indicating that radial transport of angular 
momentum by turbulent stresses may be restricted to 
interior regions near the disk midplane only |26| . 

Although MHD simulations indicate that a 
poloidal field component is essential for launching ac- 
cretion disk jets [57], it is not yet clear how such a 
field component arises in the case of accretion onto 
black holes, which have no intrinsic magnetic moment. 
The magnetic flux in the disk must originate from the 
random field in the interstellar medium (or from the 
companion star in the case of a black hole X-ray bi- 
nary). This field is then amplified into a dominant 
toroidal configuration in the disk by the MRI. Simu- 
lations have yet to demonstrate the feasibility of cre- 
ating large-scale fields via an inverse cascade process 
involving cither stochastic rcconncction of turbulent 
fields 28 or rcconncction of buoyant flux loops emerg- 
ing from the disk surface [53] ■ Models which require a 
priori large-scale flux loops [301 EH 131] , a poloidal B 



field (e.g. [IS]) or a spinning black hole ([HI [33]) to 

produce jets are too restrictive to explain the observed 
ubiquity of jets and outflows across the wide range of 
accreting systems. One of the most interesting numer- 
ical results to date is that of Machida and Matsumoto 
[19) . Their global 3D simulations show the evolution 
of a large-scale poloidal field from an initially weak 
toroidal field. However, no simulation to date has 
shown an initial random field in a generic accreting 
system evolve into a configuration favourable for jet 
production. Similarly, while Machida et al. [34] have 
studied the effect of radiative cooling on optically thin 
black hole accretion flows, no simulations to date have 
shown the effect of optically thick radiative cooling on 
the global, 3D evolution of the magnetic fields. 

We are developing new global 3D MHD simu- 
lations to directly confront this challenging problem 
and obtain new insights into the evolution of mag- 
netic field topology in black hole accretion flows. We 
aim to test our hypothesis that turbulent rcconncction 
and resistive dissipation, as well as radiative losses by 
the plasma, play pivotal roles in the evolution and 
steady-state properties of MHD accretion flows. We 
suggest that the microphysics of MHD accretion can 
govern large-scale, macrophysical phenomena that ul- 
timately determine the observational appearance and 
hence, classification, of accreting black hole sources. 

Numerical approaches to date have been limited 
by one or more of the following drawbacks: the use 
of the shearing box approximation (see, e.g., |35l I36j 
for some limitations of this approach); use of a non- 
conservative numerical scheme (e.g. [37]) resulting in 
unphysical levels of numerical dissipation which pre- 
vent a quantitative analysis of energy transport; ne- 
glect of finite resistivity making it difficult to realis- 
tically model magnetic rcconncction; and neglect of 
radiation, resulting in the unphysical situation where 
heat generated in the disk at the end of a turbulent 
cascade cannot be radiated away, so the disk puffs up. 
This also makes it impossible to compare model results 
against observations of luminosity and spectrum. 

To implement our model, we are using and ex- 
tending FLAShU [28], the public MHD code devel- 
oped at the University of Chicago. FLASH provides 
the following features which make it well suited to 
our model: it implements adaptive mesh refinement 
(AMR); it solves the equations of MHD in conserva- 
tion form thus explicitly conserving energy; it uses a 
modified piecewise-parabolic method (PPM [351 140j ) 
which is significantly more accurate than some other 
widely used codes (e.g. [3Z]); it uses the constrained 
transport method [5T] to enforce divergence free mag- 
netic fields; and it is modular and extensible, allowing 
us to develop a radiation MHD module. As such it 
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Fig. 2 A schematic illustration of the inner regions of an MHD accretion flow around a black hole. j3 is the usual plasma 
beta (the ratio of gas to magnetic pressures). The wind can drive mass loss from the disk, while the jet may be 
dominated by Poynting flux. 



represents the next generation of MHD codes. 

The organization of this paper is as follows. In 
§2, we briefly review an analytic model for turbulent 
MHD black hole disk accretion that forms the theoret- 
ical basis for our numerical investigations. In §3, we 
describe the simulations we are currently developing, 
including non-ideal MHD effects and radiation. We 
present some concluding remarks in §4. 

2. MHD Disk Accretion Theory 

The analytical foundation for our approach is 
based on the model of Kuncic and Bicknell [H] . This 
model employs a mass-weighted statistical averaging 
of the MHD equations to obtain a mean-field descrip- 
tion of turbulent MHD disk accretion that is steady- 
state and axisymmetric in the mean. Angular mo- 
mentum and energy are transported radially outwards 
by turbulent Maxwell stresses and vertically outwards 
by a large-scale MHD wind and/or jet. The inner re- 
gion of the black hole accretion disk is also surrounded 
by and magnetically coupled to a hot, diffuse corona, 
analogous to the solar corona. 

Two important observational predictions of the 
model are: 1. The disk emission spectrum is de- 
graded by the electromagnetic extraction of gravita- 
tional binding energy from the accreting matter (see 
[43l [44]); and 2. The presence of a magnetized jet 
substantially enhances the rate of mass accretion in 
the disk and hence, the rate of black hole growth, re- 
sulting in a correlation between black hole mass and 
radio emission (see [H]). The model is schematically 
illustrated in Fig. [21 The main results pertinent to our 
simulations are summarised below - more details can 
be found in the original paper [4"2"] . 



2.1 Statistical Averaging 

In the mass-weighted statistical averaging ap- 
proach, all variables are decomposed into mean and 
fluctuating parts, with intensive variables such as ve- 
locity v mass averaged according to 

Vi = Vi +v' i: (pv' t ) = (1) 

while extensive variables such as density p, pressure p 
and magnetic field B, are averaged the following way: 

p = p + p', (p')=0 , (2) 
p = p + p', (p') = , (3) 
Bi = B t + B't, (Bi) = . (4) 

For clarity, in the following equations the tilde and 
bar have been omitted from the averaged intensive and 
extensive variables, respectively: averaged quantities 
are implicitly assumed. We consider below only the 
simplest case where B» = 0, although it is straight- 
forward to generalize to the case with a nonzero net 
mean magnetic field. We have also omitted negligible 
correlation terms H 

2.2 Mass Transfer 

Integration of the mean-field continuity equation 
gives 

M a (r) + M w (r) = constant = M , (5) 

where M a (r) is the mass accretion rate and M w (r) is 
the mass outflow rate associated with a mean vertical 

2 In particular, triple correlation terms of the form (tijVj) are 
assumed negligible compared to analogous correlations with the 
mean fluid velocity (tij)v'y 



velocity at the disk surface, i.e. at the base of a disk 
wind. Under steady-state conditions, the radial mass 
inflow decreases towards small r at the same rate as 
the vertical mass outflow increases in order to main- 
tain a constant net mass flux, M, which is the net 
accretion rate at r = oo. 

2.3 Momentum Transfer 

Statistical averaging of the momentum equation 
yields 



d(pvj) d(pVjVj) 
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where 4>g = —GAI(r 2 + z 2 ) -1 / 2 is the gravitational 
potential of the central mass M, ifj = —(pv^Vj) is the 
Reynolds stress, and the turbulent Maxwell stress is 
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Integration of the azimuthal component of the 
momentum equation yields 

M a v^r - M a (n)^(n)n 

= - 2nr 2 T r( j, + 2?:r 2 T r ^(r i ) 
AM-. 



dr 



dr, (8) 



where quantities calculated at the disk surface are de- 
noted by a '+' superscript, n denotes the radius at 
the innermost stable circular orbit and 
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is the r0 component of the Maxwell stress integrated 
over the vertical scaleheight h. 

The left hand side of eqn. (jSJ) describes the change 
in angular momentum flux associated with inflow from 
an outer radius r to r\. The first two terms on the 
right hand side describe the rate of radial transport 
of angular momentum due to MHD stresses in the 
disk. The terms in the integrand on the right hand 
side describe the vertical transport of angular momen- 
tum resulting from mass loss in a wind and an MHD 
torque on the disk surface, respectively. These effects 
are not modelled in standard accretion disk theory. 
In summary, the model we consider includes contri- 
butions from both radial and vertical transport of an- 
gular momentum to the overall mass accretion rate: 
angular momentum is transported radially outwards 
by internal MHD stresses and vertically outwards by 
both mass outflows and MHD stresses acting over the 
disk surface. 



Previous simulations jT^l [T71 H5] as well as semi- 
analytic models (e.g., [25]) show that in the presence 
of a large-scale, open, mean magnetic field, angular 
momentum is transported at small disk radii more ef- 
ficiently in the vertical direction by large-scale mag- 
netic torques (Poynting flux) than radially by MHD 
turbulence. Our simulations will compare the contri- 
bution from these processes as well as from the vertical 
transport of angular momentum by mass outflows. 

2.4 Energy Transfer 

Accretion extracts gravitational binding energy 
from the accreting matter and converts it into mechan- 
ical (e.g. kinetic, Poynting flux) and non-mechanical 
(e.g. radiative) forms. The rate at which this occurs 
is determined by the keplerian shear in the bulk flow, 
s r0 = ^rdfl/dr, with dQ/dr = -§fi/r. 

The rate per unit disk surface area at which en- 
ergy is emitted in the form of clcctromagentic radia- 
tion is determined by the internal energy equation: 
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where u is the gas plus radiation energy density, F 
is the radiative flux, J is the current density, and tjj 
is the viscous stress tensor. The terms on the left 
hand side describe the total rate of change of gas plus 
radiation energy density and the terms on the right 
hand side describe work done by compression in the 
flow against the gas and radiation pressure, radiative 
losses, mean field ohmic heating, and viscous dissipa- 
tion (heating). The last term requires some comment 
since the molecular viscosity in the mean flow is gener- 
ally considered negligible in accretion disks. However, 
at the high-wavenumber end of a turbulent cascade, 
it can become important in converting the turbulent 
energy into heat. The viscous stress tensor is 



= 2vps t 



(11) 



and Sij is the shear tensor 



where v is the coefficient of kinematic shear viscosity 

(12) 
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The source terms determine the rate at which en- 
ergy is converted into random particle energy (some of 
which is then converted into radiation) and into bulk 
kinetic energy. If there are negligible changes in the 
internal energy of the gas and the turbulent energy 
is dissipated at the end of a turbulent cascade at a 
rate equivalent to its production, then eqn. (jTUJ) im- 
plies that the disk radiative flux emerging from the 



disk surface is 



-T r <f,r 



dr 
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The level of these turbulent MHD stresses avail- 
able to dissipate the internal energy in turn depends 
on how efficiently the gravitational binding energy 
extracted by accretion is converted into other forms 
(both mechanical and non-mechanical). That is, T r § 
is determined by the angular momentum conservation 
equation (c.f. |5J): 
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Substituting this into (fl~3|) yields the following more 
general solution for the disk radiative flux: 
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This is the generalized solution for the radiative 
flux of a turbulent MHD accretion disk. It can be 
expressed as 

FUr) ^ GM ^/ T) \Ur) ~ U(r)\ , (16) 



where 
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1 - 
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is a dimcnsionlcss factor that parameterizes the avail- 
able accretion energy flux, the last term describing the 
rate at which MHD stresses at the innermost stable 
circular orbit locally dissipate turbulent energy, and 
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is the fractional rate of vertical energy transport from 
the disk (the 'w' subscript denoting a wind). This is a 



correction factor which takes into account partitioning 
of accretion power into non-radiative forms. 

Note the difference between this model and stan- 
dard disk accretion theory pQ. In the latter, all the 
gravitational binding energy is locally dissipated and 
assumed to be converted to radiation: dM a /dr = 
and (t < p z ) + = so that / w = 0. This difference will be 
manifested by a disk spectrum which differs from that 
predicted by the standard model since the local disk 
temperature T(r) is reduced if energy is channelled 
away by outflows from the disk surface. This will af- 
fect the emission spectrum arising from the innermost 
regions of the disk, where the temperature is highest 
and where jets and outflows originate. Assuming local 
blackbody emission, the disk luminosity spectrum can 
be calculated by summing up the contributions from 
each annulus: 



L d , u = 2 / 7rfl„[T(r)]27rrdr, 



(19) 



where B v is the Planck function, T{r) = [F^{r)/u] 1 ' i 
is the effective disk temperature of each annulus and 
a is the Stefan-Boltzmann constant. We expect the 
disk radiative efficiency to be lower than the canonical 
~ 10% predicted by the standard model when vertical 
transport of angular momentum is important. The 
radiative efficiency is given by the ratio of disk lumi- 
nosity to accretion power: Ld/Pn- The total accretion 
power is calculated from 



^^/a^rdr. (20) 



If there is no wind mass loss from the disk, so that 
M a is constant, this reduces to the familiar result 
P a = \GMMJn YjAiaC 2 , in the Newtonian ap- 
proximation for a nonrotating black hole. 



3. MHD Accretion Simulations 

As described earlier, there is a real need for a new 
program of MHD simulations to advance our knowl- 
edge of accreting black hole systems. Our numerical 
work is motivated by the following: 

1. We need to explain the macrophysics of observed 
phenomena in AGN and other accreting systems, 
viz., high mass accretion rates and jets/winds, 
and test the hypothesis that they are related by 
large-scale MHD processes. 

2. We need to improve our understanding of the mi- 
crophysics in order to explain how small-scale, 
local MHD processes can evolve into large-scale, 
global phenomena. 

3. We need to explicitly calculate the radiation emit- 
ted by a black hole accretion disk in order to di- 
rectly compare against the observational data. 



FLASH solves the time-dependent equations of 
compressible non-ideal MHD. In non-dimensional con- 
servation form these are: 
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is the specific total energy, e is the specific internal 
energy, T v is the viscous stress tensor, g is the gravita- 
tional force per unit mass, k is the heat conductivity, 
and 77 is the resistivity. 

FLASH implements a Direct Eulcrian PPM solver 
[3l?l HO] . The constrained transport method [JT] is 
used to enforce divergence-free magnetic fields. 

3.1 Global 3D Simulations 

MHD simulations of the MRI in accretion disks 
are often conducted in a shearing box approximation 
due to the high resolution required to model MHD tur- 
bulence. However, this approach can introduce com- 
plications and pitfalls including a limited spatial scale 
for the simulations, side-effects of the shearing box 
symmetry and artifacts from the application of pe- 
riodic boundary conditions [35] as well as an aspect 
ratio dependence for MRI channel solutions [35] • We 
are now at a stage where high resolution global 3D 
simulations are possible and this is the approach we 
will take. 

3.2 Non-ideal MHD 

In the disk model of Kuncic and Bickncll [32], jets 
and/or winds are primarily responsible for transport- 
ing the angular momentum necessary for accretion to 
proceed at a rate consistent with observations of the 
most powerful astrophysical sources. We will test the 
hypothesis that the large-scale poloidal magnetic fields 
necessary for these outflows may be self-consistently 
generated in the accretion flow. Recent simulations 
[T7| show that in the presence of an externally applied 
large-scale magnetic field, angular momentum trans- 
port by the vertical (<pz) Maxwell stress is compara- 



ble to its radial (rip) component. Magnetic reconnec- 
tions can have a significant influence on magnetic field 
topologies [37] US] . Notwithstanding the high degree 
of ionisation of plasmas in the accretion disks of AGN 
and X-ray binaries, we suggest that these reconnec- 
tions and non-ideal MHD effects in general cannot be 
neglected. Even in numerical models that do not ex- 
plicitly include non-ideal effects, they can appear in 
the form of numerical resistivity which is difficult to 
control and quantify. By explicitly modelling a finite 
resistivity, we will explore its effect on the evolution 
of the magnetic field topology, particularly the emer- 
gence of a significant z-component which is necessary 
to produce high mass accretion rates. 



3.3 Radiation 

The inclusion of radiation in our simulations is 
imperative for directly comparing to the observational 
data, which is almost exclusively in the form of pho- 
tons detected in various wavebands. Most of the emis- 
sion that characterises quasars and other AGN is at- 
tributed to the putative accretion disk and peaks at 
optical-ultraviolet spectral energies. Radiative trans- 
fer will also be required to transport the internal en- 
ergy dissipated in the disk plasma at the end of a 
turbulent cascade, i.e. to cool the disk. To date, no 
simulations have investigated the effect of optically 
thick radiative cooling on the MRI in full global 3D. 
Including radiation will also allow us to test whether 
the blackbody emission from the disk is modified by 
outflows. In addition, it may be that regions of the 
disk where radiation dominates may be thermally un- 
stable [12], thus affecting the dynamics. 

Implementing radiation is computationally very 
demanding; numerically solving the full radiative 
transfer problem in 3D is currently not feasible. In- 
stead, a common approach is to average over frequency 
and solve the equations in the flux-limited diffusion 
(FLD) approximation [50] [5T], a technique which still 
allows one to approximate the emergent spectrum. A 
"flux hmiter" is used to interpolate between the opti- 
cally thin and optically thick cases, giving a reasonable 
measure of the energy carried away by radiation [52] . 
FLD has previously been implemented in a shearing 
box in a reference frame co-moving with the fluid (e.g. 
53, 37J). Simulations show a stratified disk in contrast 
to the standard model [54], and that radiative diffu- 
sion dominates Poynting flux throughout the disk and 
the upper layers are magnetically supported and inho- 
mogeneous, likely affecting the emergent thermal spec- 
trum [55J [5S [57] . The implementation of MHD with 
FLD described in [55] achieves energy conservation by 
using a mixed-frame numerical scheme, evaluating ra- 
diation quantities in the lab frame and fluid opacities 
in the co-moving frame, thus enabling us to address 



radiation in a quantitative way. The algorithm also 
provides improved speed and accuracy compared to 

[531 S3- 

Our goal is to calculate the steady-state emission 
spectrum of a turbulent MHD disk around a super- 
massive black hole in the nucleus of a galaxy and to 
compare the predicted spectrum with the observed op- 
tical spectra of quasars (see, e.g., [43]). 

4. Concluding Remarks 

The publication of these proceedings coincides 
with the 50th anniversary of the discovery of the MRI 
[3], which has had such a profound impact on our 
understanding not only of accretion disks - arguably 
nature's most powerful energy source - but plasmas 
in general, ranging from laboratory scales to galaxy 
scales. Further landmarks in accretion disk theory 
came with the laying down of the standard theory [T] 
in what remains the most cited paper in all of astro- 
physics; the discovery that large scale magnetic fields 
can vertically transport matter, energy and momen- 
tum [24] ; the first numerical MHD simulations [9] [10] ; 
as well as the rediscovery of the MRI in accretion disks 
and the demonstration by numerical simulations that 
this is indeed the MHD process anticipated by the an- 
alytic standard theory necessary to produce the turbu- 
lent radial transport of energy and momentum [SJ [5] . 
Despite these major steps forward in assembling the 
components of a complete accretion disk theory, we 
have not to date seen numerical simulations which 
can self- consistently produce all the salient features of 
quasars and other high energy astrophysical sources: 
high mass accretion rates, outflows, winds and jets, 
the formation of a magnetised corona, and the ob- 
served thermal spectrum. 

We expect it will again be magnetic fields which 
will hold the key to resolving these outstanding is- 
sues. To this end, we are developing new 3D global 
MHD simulations to test the hypothesis that small- 
scale, stochastic fields can self-consistently generate 
the large-scale poloidal magnetic fields necessary for 
the transport of energy and momentum from accret- 
ing matter necessary to produce each of the above 
observed features. The analytical basis for this gener- 
alized model was laid down by the companion work 
|42] to our current numerical modelling, which the 
rapid advances in computing power accompanied by 
new codes and more efficient algorithms has now made 
possible. 

It is an exciting time for accretion disk theory and 
for our understanding of the microphysics that drives 
these and other plasma systems in nature. 
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